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ABSTRACT 


This  report  aids  in  the  calculation  of  blackbody  radiation.  The 
equations  of  blackbody  radiation  and  photon  emission  are  summarized 
and  a  computer  program  is  described  that  is  useful  for  calculations 
concerning  blackbodies.  The  input  to  the  program  consists  of  the 
units  desired,  the  temperature,  and  wavelengths.  The  output  includes 
the  physical  constants  in  the  desired  units  and  the  evaluation  of  the 
equations  for  the  temperature  and  wavelengths  specified.  The  equa¬ 
tions  for  photon  and  power  emission  are  written  in  terms  of  Debye 
functions.  The  computer  program  uses  a  subroutine  that  calculates 
to  five  significant  figures  the  Debye  functions  of  orders  from  one 
through  six.  The  method  and  error  analysis  of  this  calculation  are 
included. 
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NOMENCLATURE 

Conversion  factor  (Eq.  (10)) 

Term  of  infinite  series  (Eq.  (49)) 

Ratio  of  over  Ta 
Ratio  of  i/Me  over  Ta 
Product  of  and  Ta 

Product  of  and  Ta 

Bernoulli  numbers 
Speed  of  light 
First  radiation  constant 
Second  radiation  constant 
Debye  function 

Debye  function  or  its  complement,  depending  on  i 

Monochromatic  emissive  power  in  terms  of  wavelength 

Monochromatic  emissive  power  in  terms  of  frequency 

Total  emissive  power 

Emissive  power  in  terms  of  wavelength 

Emissive  power  in  terms  of  frequency 

Relative  error 

Monochromatic  photon  emittance  in  terms  of  wavelength 

Monochromatic  photon  emittance  in  terms  of  frequency 

Total  photon  emittance 
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Photon  emittance  in  terms  of  frequency 

Planck  constant 

Boltzmann  constant 

Subscript  indexing  terms  in  infinite  series 
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j?  Superscript  distinguishing  different  series 

m  Number  of  terms  of  infinite  series  retained 

n  Order  of  Debye  function 

r  Defined  in  Eq.  (57) 

£ 

Rnm  Remainder 

Sn  Infinite  series 

a 

Snm  Partial  sum 

T  Temperature 

Ta  Absolute  temperature 

x  Defined  in  Eq.  (20)  or  independent  variable  of  Debye  function 

£  Zeta  function 

X  Wavelength 

^Me  Wavelength  of  maximum  e^ 

Xj/[f  Wavelength  of  maximum  f^ 

v  Frequency 

^ IVLe  Frequency  of  maximum  e^ 

Frequency  of  maximum  fbu 
o  Stefan-Boltzmann  constant 

t  Ratio  of  Ffo  to  Ta3 
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SECTION  I  . 
INTRODUCTION 


The  blackbody  concept  is  important  in  theoretical  and  applied 
physics  and  in  radiation  heat  transfer.  Some  of  its  properties,  such  as 
the  Stefan-Boltzmann  law  and  Weins  displacement  law  can  be  derived 
from  the  second  law  of  thermodynamics.  The  search  for  a  theoretical 
derivation  of  the  spectral  distribution  function  of  blackbody  radiation, 
which  fits  experimental  data,  led  to  Planck's  Radiation  law,  which 
was  the  birth  of  quantum  mechanics  (Ref.  1).  The  complete  knowledge 
of  the  spectral  distribution  of  blackbody  radiation  makes  it  a  valuable 
standard  for  checking  the  performance  of  wavelength-dependent  devices 
such  as  optical  filters,  spectrometers,  and  radiometers.  Many  calcu¬ 
lations  of  radiation  heat  transfer  are  based  on  the  assumption  that  the 
surfaces  are  black  or  gray. 

In  certain  areas,  calculations  involving  blackbody  radiation  are 
frequently  required.  These  calculations  necessitate  looking  up  or 
remembering  the  equations  and  physical  constants,  performing  unit 
conversions,  using  tables,  and  evaluating  the  equations.  There  are 
several  useful  aids  in  the  calculations,  such  as  the  photon  slide  rule, 
the  radiation  slide  rule,  and  tables  of  Refs.  2,  ,3,  and  4,  respectively. 
This  report  and  the  described  computer  program  should  also  be  useful. 
A  summary  of  the  equations  of  blackbody  radiation  and  photon  emission 
is  included.  The  input  to  the  computer  program  includes  the  units  de¬ 
sired,  the  temperature,  and  the  wavelengths;  the  output  consists  of  the 
physical  constants  in  the  desired  units  and  the  evaluation  of  the  equa¬ 
tions  for  the  temperature  and  wavelengths  specified. 

The  photon  and  power  emissions  between  wavelengths  are  written 
in  terms  of  Debye  functions.  The  computer  program  uses  a  subroutine 
that  calculates  to  five  significant  figures  the  Debye  functions  of  orders 
from  one  through  six.  The  method  and  error  analysis  of  this  calcula¬ 
tion  are  good  examples  of  numerical  application  of  infinite  series;  this 
analysis  is  given  in  Section  IV,  independent  of  the  other  material. 


SECTION  II 
EQUATIONS 


The  energy  emitted  by  a  blackbody  per  unit  time  and  per  unit  area 
in  a  frequency  range  du  is  e^  du  where  e^  is  the  monochromatic 
emissive  power  in  terms  of  frequency  and  is  given  by 
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ebt/ 


2n-  h  v 3 


(1) 


The  energy  emitted  per  unit  time  and  per  unit  area  in  a  frequency  range 
from  0  to  v  is  the  emissive  power  in  terms  of  frequency  and  is 


v 

Eby  ~  ^  ebv  dv 


(2) 


The  total  energy  emitted  per  unit  time  and  per  unit  area  from  all  fre¬ 
quencies  is  the  total  emissive  power 


eby  di/ 


(3) 


The  maximum  value  of  eby  occurs  at  the  frequency  ^jy[e  which  is  the 
nonzero  root  of 


dv 


0 


(4) 


A  photon  is  a  packet  of  energy,  h v.  Thus,  the  number  of  photons 
emitted  by  a  blackbody  per  unit  time  and  per  unit  area  in  a  frequency 
range  du  is 

hv  Al/  =  Av  (5) 


and  the  monochromatic  photon  emittance  in  terms  of  frequency,  fbi,,  is 


or 


fbv  = 


CbV 

hv 


(6a) 


f  2  n  v2 _ 

clexpG^7)  -  !]  (6b) 

The  number  of  photons  emitted  per  unit  time  and  per  unit  area  in  a 
frequency  range  from  0  to  v  is  the  photon  emittance  in  terms  of  fre¬ 
quency  and  is 

Fbv  =  J  fhv  di/  . 


The  total  number  of  photons  emitted  per  unit  time  and  per  unit  area 
from  all  frequencies  is  the  total  photon  emittance  and  is 


(8) 
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The  maximum  value  of  ft,*/  occurs  at  the  frequency  which  is  the 
nonzero  root  of 


IhiL.  =  o 

dv 


(9) 


Wavelength  and  frequency  are  related  by 


(10) 


where  a  is  a  conversion  factor,  such  as  cm  per  At,  since  sometimes 
different  length  units  are  used  for  c  and  X.  The  energy  emitted  by  a 
blackbody  per  unit  time  and  per  unit  area  in  a  wavelength  range  dX  is 
ebX  d''-»  where  e^^  is  the  monochromatic  emissive  power  in  terms  of 
wavelength  and  is  given  by 


eb\  =  -ebv 


dv 

dA 


or 


ebX  = 


e4A*  exp  ( - ^  -  1 

L  \akAT  B  /  J 


which  is  commonly  written  as 


ebX  = 


where 


and 


c,  = 


2rrhc2 


he 

ak 


(11a) 

(lib) 

(He) 

(12) 

(13) 


The  energy  emitted  per  unit  time  and  per  unit  area  in  a  wavelength 
range  from  0  to  X  is  the  emissive  power  in  terms  of  wavelength  and  is 

EbX  =  J  ebXdX  (14) 

o 


The  maximum  value  of  e^  occurs  at  the  wavelength  X]y[e  which  is  the 
nonzero  root  of 


debA 


=  0 


(15) 


Note  that  Xjyje  is  not  the  wavelength  corresponding  to  ^]vie. 
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The  number  of  photons  emitted  per  unit  time  and  per  unit  area  in  a 
wavelength  range  dX  is  fbX^*  where  fbA  is  the  monochromatic  photon 
emittance  in  terms  of  wavelength  and  is 

fbA  =  -fb„  ~  ( 16a) 

which  may  be  written  as 


fbA 


2nc 


hfe)- : 


(16b) 


The  number  of  photons  emitted  per  unit  time  and  per  unit  area  in  a 
wavelength  range  from  0  to  X  is  the  photon  emittance  in  terms  of  wave¬ 
length  and  is 


fbA  <fX 


(17) 


The  maximum  value  of  f^  occurs  at  the  wavelength  Xjy[f  which  is  the 
nonzero  root  of 


dfbA 

dA 


=  0 


(18) 


Note  that  Xjy[f  is  not  the  wavelength  corresponding  to  rjy[f. 


The  emissive  power  and  photon  emittance  can  be  written  in  terms 
of  Debye  functions  (Ref.  5).  A  table  of  Debye  functions  is  given  in 
Ref.  5,  and  a  method  to  calculate  them  is  given  in  Section  IV  of  this 
report.  The  nth  order  Debye  function  is  defined  as 


-  r^br-  * 


making  the  substitution 


hi/ 

kT„ 


Combination  of  Eqs.  (1),  (2),  and  (19)  results  in 

Eb„  =  T a  DaW 

c  h 

The  substitution  into  Eq.  (7)  results  in 


Fbv  =  Ta3  D2(x) 

c  a 


(19) 


(20) 


(21) 


(22) 
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The  value  of  the  Debye  functions  at  infinity  (Ref.  4)  is 


D„(~)  =  n!  tfn  +  1) 

(23) 

where  the  zeta  function  is  defined  as 

<(n>  -  £  T 5r 

k=l 

(24) 

Using  this  notation,  and  Eqs.  (3)  and  (21),  one  obtains 

Eb  =  <(4)  Ta‘ 

c  h 

(25a) 

which  is  commonly  written  as 

Eb  =  oTaA 

(25b) 

The  value  of  £(4)  is  (Ref.  5) 

£(4)  =  JL.  =  1.08232  32337  11138  19152 

(26) 

thus,  cr  is 

<x'=  2*V 

15cah3 

(27) 

From  Eqs.  (8)  and  (22)  one  obtains 

Fb  =  -?■£-  £(3>  Ta3 

cJ  h3 

(28) 

which  can  be  written  as 

Fh  =  rTa3 

(29) 

where 

'  =  %  . 

(30) 

There  is  no  exact  expression  for  £(3)  as  there  was  for  ?(4); 
its  value  is  (Ref.  5) 

however. 

C(3)  =  1.20205  69031  59594  28540 

(31) 

The  substitution  of  Eq.  (20)  into  Eq.  (14)  results  in 

FbA  =  Eb  —  Ebi/ 

(32) 

and  similarly 

FbA  =  Fb  -  Fbi/ 

(33) 
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Performing  the  differentiation  of  Eqs.  (4),  (9),  (15),  and  (18),  the 
result  in  each  case  is  of  the  form 


X  =  n(l  -  e-*)  (34) 

The  nonzero  roots,  xn,  of  Eq.  (34)  are  given  in  Table  I  for  the  values 
of  n  from  2  through  10.  Thus  is  obtained 


'•'Me  = 

vm  = 
Am  e  T  a  = 
AMfTa  = 

where 

b,  = 
b2  = 
b,  = 
b4  = 


b,  Ta 

(35) 

b2  Ta 

(36) 

bs 

(37) 

b« 

(38) 

k  x3 

h 

(39) 

k  x2 

h 

(40) 

he 

akxs 

(41) 

he 

akx. 

(42) 

The  equations  of  this  section  summarize  the  expression  for  the 
quantities  of  blackbody  radiation  and  photon  emission  which  are  most 
frequently  to  be  calculated.  The  next  section  describes  a  computer 
program  that  relieves  the  calculator  of  the  labor  of  finding  physical 
constants  and  conversion  factors,  performing  unit  conversions,  using 
tables,  and  evaluating  the  equations. 


TABLE  I 

ROOTS  OF  x  =  n  (1  -e-*) 


n 

xn 

2 

1. 59362  42601 

3 

2. 82143  93721 

4 

3. 92069  03949 

5 

4. 96511  42317 

6 

5. 98490  12264 

7 

6. 99357  56867 

8 

7.  99730  90676 

9 

8. 99888  80761 

10 

9.  99954  57944 
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SECTION  III 
COMPUTER  PROGRAM 


The  computer  program  was  written  in  IBM  System/ 360  FORTRAN  IV 
language.  Its  calculations  can  be  divided  into  three  groups: 

1.  Unit  conversions  and  calculation  of  the  physical  constants. 

2.  Calculation  of  quantities  dependent  on  temperature,  but  not 
wavelength. 

3.  Calculation  of  wavelength -dependent  quantities. 

The  flow  chart  of  the  program  is  given  in  Fig.  1,  and  a  description  of 
the  logic  is  given  below  followed  by  an  explanation  of  the  input  and  output. 


Fig.  1  Flow  Chart  of  the  Computer  Program 
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Fig.  1  Continued 
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Fig.  1  Continued 
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Fig.  1  Continued 
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Calculate  Wavelength-Dependent  Quantities 


efev  -  2  it  h  v3/[c2(exp  x  -  1)  ] 
fbv  -  2  T  v2/[c2 (exp  x  -  1)] 
ebx  =  2  t  h  c2/[a4  X5(exp  x  -  1)  ] 
fbx  =  2  t  c/[a3  X4(exp  x  -  1) ] 

‘  2  '  *"  D3U)/1q2  h3> 

rbv  -  2  T  “3  i  D2(*)/<=2  »3> 

Eb\  Eb  ®bv 


Fig.  1  Concluded 
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One  good  feature  of  the  program  is  the  simplicity  of  its  logic.  In¬ 
put  is  via  one  READ  statement: 

READ  IW,  IL,  IS,  IT,  IE,  L,  N,  T,  QO,  DQ 

According  to  the  format: 

616, 

16,  3E12.  0 

The  input  variables  are: 

IW 

Indicator  of  the  length  unit  of  X 

IL 

Indicator  of  the  length  unit  of  c 

IS 

Indicator  of  the  time  unit 

IT 

Indicator  of  the  temperature  unit 

IE 

Indicator  of  the  energy  unit 

L 

Indicator  of  whether  X,  v ,  or  x  is  the  independent 
variable 

N 

Number  of  values  of  the  independent  variable  for  which 
the  wavelength  dependent  quantities  are  to  be  calculated 

T 

Temperature 

QO 

Starting  value  of  the  independent  variable 

DQ 

Increment  of  the  independent  variable 

The  program  reads  a  card  and  performs  the  indicated  computations  ex¬ 
plained  below.  That  is  one  complete  cycle  of  the  program.  The  program 
then  branches  to  the  start  and  begins  another  cycle  by  reading  another 
card.  Execution  is  terminated  when  the  end  of  the  data  is  reached. 

The  first  calculations  are  the  unit  conversions  and  the  calculation 
of  the  physical  constants.  The  basic  physical  constants  are  a,  c,  h,  k, 
and  Tref.  All  the  other  constants,  a,  t,  b2,  b3,  b4,  b5,  cj,  and  C2, 
are  expressed  in  terms  of  the  basic  constants.  The  program  contains 
the  values  of  the  basic  constants  in  the  International  System  of  Units  (SI) 
(Ref.  6)  and  the  values  of  various  conversion  factors  (Table  II).  Unit 
conversions  indicated  by  the  first  five  input  variables  are  performed  on 
the  four  basic  constants,  then  the  other  constants  are  calculated  by  their 
defining  equations.  The  valid  values  for  the  first  five  input  variables 
and  the  units  which  the  values  indicate  are  given  in  Table  II.  This  table 
should  be  referred  to  when  preparing  input  data  for  the  program.  If  IW 
is.  left  blank  on  an  input  card,  then  this  first  group  of  calculations  is 
omitted.  The  values  of  the  physical  constants  used  in  the  previous  cycle 
are  retained  for  the  second  and  third  group  of  calculations. 
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TABLE  II 

VALUES  OF  BASIC  PHYSICAL  CONSTANTS  AND  CONVERSION  FACTORS 


Values  of  Basic  Physical  Constants 
in  SI  Units  (Ref.  5) 


a  =  1  m/m 

c  =  2.997925  x  108  m/s 
h  =  6.6256  x  1(T34  Js 
k  =  1.  38054  x  10*23  J/°K 


IVV  and  IL 


Length  Units  (Ref.  5) 


Units  Indicated 

Conversion  Factor 

/J 

1  X  10'6  m/a 

cm 

0.  01  m/cm 

m 

1  m/m 

in. 

0.  0254  m/in. 

ft 

0.  3048  m/ft 

Time  Units  (Ref. 

5) 

IT 

Units  Indicated 

Conversion  Factor 

i 

s 

1  s/s 

2 

min 

60  s/min 

3 

hr 

3600  s/hr 

Temperature  Units  (Ref.  5) 

IT 

Units  Indicated 

Conversion  Factor 

1 

°K 

1  'K/’K,  T ref  -  0 

2 

•R 

0.5555556  °K/°R,  Tref  =  0 

3 

°c 

1  “K/'C,  Tref  =  273.  15 

.4 

“F 

0.5555556  'K/T,  Tref  =  459.67  ] 

Energy  Units  (Ref.  5) 

IE 

Units  Indicated 

Conversion  Factor 

1 

erg 

1  x  10"7  J/erg 

2 

J 

1  J/J 

3 

cal  (mean) 

4.  19002  J/cal 

4 

ft -lb 

1.  35582  J /ft -lb 

5 

Btu  (mean) 

1055.  87  J/Btu 
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The  quantities  Ta,  Eb,  Fb,  ^Mf»  %e>  and  ^-Mf  depend  on  the 

temperature,  but  not  upon  the  wavelength.  These  are  the  quantities 
calculated  by  the  second  group  of  calculations.  If  T  is  less  than  or  equal 
to  absolute  zero,  then  the  second  group  of  calculations  is  omitted.  The 
values  of  the  temperature -dependent  quantities  used  in  the  previous 
cycle  are  retained  for  use  in  the  third  group  of  calculations. 

The  wavelength -dependent  parameters  are  X ,  u,  x,  e^v,  fb„ ,  ebx* 
fbX»  Ebi/»  Ebi/»  EbX*  and  Eb^.  These  are  the  quantities  calculated  by 
the  third  group  of  calculations.  The  variables  X,  v ,  and  x  are  inter¬ 
dependent  and  any  one  can  be  chosen  as  the  independent  variable.  This 
choice  is  indicated  by  L,  with  values  of  1,  2,  and  3  indicating  X,  v , 
and  x,  respectively.  The  use  of  x  as  the  independent  variable  has  the 
advantage  that  the  principal  part  of  the  spectrum  is  in  the  range 
0.  1  <  x  <  10,  whereas  the  corresponding  range  for  X  and  v  depends  on 
the  temperature.  The  wavelength-dependent  variables  are  calculated 
for  N  values  of  the  independent  variable,  starting  with  QO  and  incre¬ 
menting  by  DQ  each  time.  The  value  of  Q  is  limited  to  values  corre¬ 
sponding  to  x  greater  than  0.  001  and  less  than  150.  If  Q  is  outside  this 
range,  it  is  set  equal  to  the  nearest  limiting  value.  If  N  is  left  blank  on 
an  input  card,  then  the  third  group  of  calculations  is  omitted. 

Figure  2  illustrates  an  input  form  designed  for  this  program  and 
data  for  a  sample  run.  The  corresponding  output  is  shown  in  Fig.  3. 

The  FORTRAN  listing  of  the  computer  program  is  given  in  Appendix  I. 

A  list  of  FORTRAN  variables  is  given  in  Table  III. 


Columns 

1 

2 

3 

1 

5 

6 

1 

I 

ii 

12 

37  -  48 

N.  Input 

^\Var  lable 

Card 

No. 

IV 

IL 

IS 

IT 

IE 

0 

N 

T 

QO 

DQ 

1 

■ 

B 

3 

4 

5 

■ 

■ 

■ 

■ 

B 

B 

-500. 

2 

1 

4 

2 

2 

3 

3 

3 

3 

1 

1 

2 

__ 

4 

2 

2 

1 

3 

1 

-500. 

5 

1 

2 

1 

1 

2 

20. 

6 

77. 

7 

300. 

8 

3 

■ 

■ 

i 

1 

6000. 

0. 

1. 

9 

3 

■ 

8 

20. 

20. 

_ 

■ 

■ 

■ 

■ 

Fig.  2  Input  Data  for  Sample  Run 
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BLACK  BCD Y  EMISSION 


CJ1 


1  A 

IL 

IS 

I  I 

lb 

A 

C 

t 

K 

TREF 

TAM 

B2 

S3 

H4 

BS 

Cl 

C2 

MICK 

FI 

HR 

FOEv. 

BTU 

3.28C8E-06 

3.S4C9E 

12 

1.7431F- 

4C 

7.2633E- 

27 

4.5967E  02 

1.710Ht-O9 

&  .  7  I  9  ?c 

16 

6.64  lit 

13 

1.1753E  14 

b.  605  5E 

C3 

5.2160E  03 

1.  1851  E 

08 

2.SB9HE 

C4 

MICK 

IK 

MIN 

ROEG 

CAL 

3.S37CE-CS 

7.0817F 

11 

2.6355C- 

36 

1.82C5E- 

24 

C.O 

4.99J7E— 1 L 

1.0092c 

13 

1.1068F 

12 

r.ossoE  12 

6.6055E 

03 

5.2160E  03 

3.4566E 

06 

2.5B98F 

C  4 

4 

M 

S 

KOEG 

J 

l.OCCOb  00 

2.S579E 

C9 

6.6256E- 

34 

1.38C5E- 

23 

C.C 

5.6b*8C-0H 

1.5204E 

15 

3.3205F 

10 

5 .d  7  EOF  10 

3.O697F-03 

2.  B978F-C  3 

3.7415E- 

16 

1.4339E- 

•C? 

CM 

CM 

S 

COEG 

ERG 

l.OCCOE  00 

2.99  79F 

LC 

6.6256F- 

■27 

1 .3805F- 

16 

2.7315E  02 

5*66  jJC—05 

1.52C4E 

11 

3. J205F 

1C 

5 .3  7B9F  10 

3.6697E-01 

2.89  78E-01 

3.74I5E- 

05 

1.4388E 

CC 

MICK 

CM 

S 

KOEG 

J 

1.0CC0E-C4 

2.9979E 

10 

6 .  6256E- 

■34 

1.38C5E- 

23 

C.O 

5.66*dF“i2 

1.5204b 

11 

3.  3205F 

1C 

5.378SE  10 

3.6697E 

03 

2.8978E  03 

3.7415E 

04 

1.4383E 

C4 

I  ' 

r  a 

EB 

fb 

NME 

NMF 

LME 

IMF 

2.0000b  0L 

2.oooor 

01 

9.C716F- 

■C  7 

1.2163E  IS 

1. 1758E 

12 

6. 64 1  IE  11 

I.44H9E 

02 

1.8349E 

C2 

7.7900C  HI 

7.7000E 

01 

1.9931b- 

■04 

6.9412F  lo 

4.0267E 

12 

2.556BC  12 

3*  7634 E 

01 

4.7655E 

Cl 

3.G0C9f  02 

3.0000E 

02 

4 .5925C- 

■02 

4.1GS1E  IB 

1.7637F 

13 

9.S616E  12 

5.S5S3E 

CO 

1.2232E 

Cl 

6-CCCCL-  03 

6.9000E 

03 

7.3430E 

03 

3.2B41E  22 

J.5273E 

14 

1 .99  23t  14 

4. B257E- 

-01 

6.I162F- 

•01 

LAM.jOA 

N  L 

X 

EM 

£n 

EL 

FL 

EHN 

FRN 

EBL 

FBL 

2.39aOC  0 3 

L • 25C2L 

11 

i.ccooe- 

■0  3 

9 .04720-  IB 

1.0022E 

OS 

4.7167E-10 

5.6943E 

12 

7.  77C3E- 

■C  7 

6 . 8  2  79  E 

15 

7.3480E  03 

3.264  IE 

22 

2. 39-JOE  03 

1.2502C 

14 

i.orcot' 

ro 

3.2674C-12 

6.3591E 

07 

2. 746  It  02 

3.3153F 

21 

2 . 5437T 

C2 

4.8349E 

21 

7.0936F  03 

2.6006E 

22 

1.199CC  00 

2.5004C 

14 

2.CCCCF 

CO 

1.1  333E-11 

6 . 84C9E 

07 

2.3624E  C.3 

1.4266E 

22 

1.231  IE 

02 

1.3471E 

22 

6.0170F  03 

1 .9  3  7CE 

22 

7.09 JHE-01 

3 .75061 

14 

i.cooni 

CC 

1.23C4I-11 

5.1526E 

07 

6.0C78C  03 

2.4L77E 

22 

2.E879E 

C? 

2. 1061E 

2? 

4.46C1F  03 

1 .1780E 

22 

5*9*)  ■jOF*  -  Cl 

■5.TOC7E 

14 

4.COO0E 

f  0 

1.0BC7F-11 

3.261  RE 

C  7 

9,01 SOF  C  3 

2.7209F 

22 

4.387Cr 

03 

2.6288E 

22 

2 .96 1  IF  C 3 

6. 5527E 

21 

4, 7960L-G1 

5.2509E 

14 

5.U  DOE 

03 

7.6747E-12 

1.35  31F 

07 

1.00C3E  04 

2.4152E 

2? 

5.5443E 

C  3 

2.9426E 

22 

1.8037E  03 

3.4 151 E 

21 

3 • 99b6fc-P 1 

7.501  It 

14 

6.0C OOF 

00 

4 • 35  791"—  12 

9.7747E 

06 

9.1177E  03 

1.8346F 

22 

6. 22C5E 

Cl 

3. 1 1 46  E 

22 

1.0275E  03 

1 ■ 6  94  8E 

21 

3.425  7E-01 

•1.751  3t 

14 

7.C000E 

CO 

2.H335E-12 

4.8867E 

06 

7. 2384t  03 

1.2484E 

2? 

6.7927E 

C3 

3.2C31E 

22 

5.5532E  C2 

8. 10CCE 

ZC 

2«  99  75E-C 1 

1.000  1C 

15 

H.OrOQE 

ro 

1.5551E-12 

2. J467F 

06 

5.lfle7E  C3 

7.83C1F 

21 

7.C6C3F 

03 

3.2465E 

22 

2.8776F  C2 

3.75E2E 

20 

2.66'*4E-0l 

1. 1252E 

15 

s.onot 

CO 

E.I4  SGl-13 

1.0924E 

06 

3.4190E  03 

4.61  SIC 

21 

7.203SE 

03 

3.2671E 

22 

1 .44 1 2E-  C2 

1 .7028E 

2C 

2. 394CG-0L 

1.2502b 

15 

i.ccoof 

01 

4.10931-13 

4.9609F 

OS 

2.  1424C  C  3 

2.5864E 

21 

7.2778E 

C3 

3 .2765 E 

22 

7.0174E  Cl 

7. 5664E 

19 

1. 19 -JOE-01 

2. 50C4E 

15 

2.Cf  03F 

01 

1 .VJ24i"-I6 

9.CC86E 

01 

3. 1 1 23 F  CO 

1.8787E 

18 

7.  34BCF 

C  3 

3 • 2841 E 

22 

2.I7S0E-C2 

1.2445E 

1 6 

5  .9930E-02 

5.90C7E 

15 

4.0CC0E 

01 

2.4609E-24 

1.42  73E- 

-07 

2.0S28E-07 

6.19SSE 

10 

7. 749CE 

C3 

3. 2841 F 

22 

3 .3 19  It-  10 

9. 76 1 3E 

C  7 

3.9966F-02 

7.501  IF 

15 

6.CC00C 

or 

1.7119F-32 

3.4445E  ■ 

-15 

2.2I29C-15 

6.4648E 

02 

7.3430F 

C  3 

3.2341E 

22 

2.25C9E-18 

4.45?Zt-gi 

2.99  75C-02 

l  ,')f)0  lb 

16 

s.rrooF 

01 

8 ■ 3637F-4 1 

1.2621E 

-23 

2. 79C7T-2 3 

4.21L3E- 

-C6 

7.14B0C 

C  3 

3.2841E 

22 

1 .0858E-26 

1.6179E-09 

2.3940E-02 

1 .2502E 

16 

i.0(  oor 

02 

3.  S670E-49 

4.3648E 

-32 

1.7554E-31 

2.119ZF- 

-14 

7. 3480E 

C3 

3 .2841 E 

22 

4.3  382E-35 

5. 1  b44t- 

-10 

1.998 JE-02 

l  .59021: 

lo 

1.2C00C 

02 

l  .1992C-57 

1.206SE 

-4  C 

9.00 ;0r-40 

9.C574E-23 

7.743CF 

Cl 

3. 2341  E 

22 

1  .5372E-43 

1 . 7  I2RC-02 

l. 75C3fc 

1  b 

1.4000E 

02 

3.9250F-66 

3.3847F 

-49 

4.0108F-4H 

3.4586E- 

-31 

7. 24eCE 

C3 

3.2841E 

22 

5  .0 1 37E-52 

4.2 S23E-35 

1.593 7E-02 

1.M753E 

16 

l.srcoE 

02 

2.101 7*r —  70 

1.  7640E 

-S3 

Z.571CF-52 

2 .0692b—  35 

7. 348CE 

C  3 

3.2841E 

22 

2.7956E— 56 

2.234*tf 

Tbt  END 


Fig.  3  Output  of  Sample  Run 
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TABLE  III 

FORTRAN  VARIABLES 


FORTRAN  Variable 

Algebraic  Variade 

A 

a 

£2 

b2 

B3 

k>3 

B4 

b4 

B5 

b5 

C 

c 

CR 

c  in  m/s 

Cl 

Cl 

C2 

c2 

EB 

Eb 

EBL 

EbX 

EBN 

Ebi / 

EL 

ebX 

EN 

ebt/ 

FB 

Eb 

FBL 

EbX 

FBN 

Ebi/ 

FL 

fbx 

FN 

^b  v 

H 

h 

HR 

h  in  Js 

K 

k 

K0H 

k/h 

KR 

k  in  J/°K 

LAM 

X 

LME 

xMe 

LMF 

xMf 

NME 

‘'Me 

NMF 

"Mf 

NU 

V 

PI 

ir 

SIG 

a 

T 

T 

TA 

Ta 

TAU 

T 

X 

X 

X2 

x2 

X3 

x3 

X4 

x4  1 

X5 

x5 

ZETA3 

5(  3) 

ZETA4 

5(4) 
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SECTION  IV 

DEBYE  FUNCTION  CALCULATION 


The  Debye  function  can  be  calculated  by  truncated  infinite  series. 
Theorems  of  infinite  series  were  used  to  perform  an  error  analysis  on 
this  calculation.  Of  special  significance,  a  criterion  was  obtained  for 
the  error  that  was  independent  of  the  order  of  the  function.  This 
criterion  was  used  to  deduce  conditions  for  the  calculation  of  the  Debye 
function  to  five  significant  figures. 


4.1  DEFINITIONS  AND  NOTATION 


Adding  superscripts,  1  and  2,  to  distinguish,  respectively,  the 
Debye  function  and  its  complement,  the  Debye  function  is 

D“(x)  =  /  ;r:'T  dt  (43) 

0 

and  its  complement  is 

D»(x)  =  J  et  -  i  dt  (44) 


The  sum  of  the  Debye  function  and  its  complement  is 
Dln(x)  +  DJ,(x)  =  f  - —  dt  =  n!  <(n  +  1) 

where  the  zeta  function  is  defined  by  Eq.  (24). 

An  infinite  series  expansion  for  the  Debye  function  is 

Diw  =  xn  [—  -  _i —  +  i  B“*2k  1 

U  2(n  +  ])  "1  (2k  t  n)  (2k)! J 

for  |  x  |  <  27r  and  for  n  >  1. 


(45) 


(46) 


An  infinite  series  expansion  for  the  complement  of  the  Debye  func¬ 
tion  is 


oo  n 


DJ„(x)  =  £  £ 


n !  x1*-!  e-tx 


(n  -  j)I  k3"r  1 


(47) 


for  x  >  0  and  for  n  >  1. 
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The  B2k  Of  Eq.  (46)  are  Bernoulli  numbers  and  are  discussed  in 
Ref.  5.  Of  basic  importance  to  this  analysis  is  the  inequality 


2( 2k) ! 
<2»r)2k 


<  F2k 


_  2{2k)! 

(2ff)2k 


[ 


I 

1  -  2  1—  2k 


] 


(48) 


The  Debye  function  can  be  calculated  by  using  a  partial  sum  of  one  or 
the  other  of  the  two  series,  depending  on  the  value  of  x.  The  remaining 
function  can  be  calculated  from  Eq.  (45). 


In  the  following  analysis,  the  theorems  used  are  proved  and  dis¬ 
cussed  in  Ref.  7.  The  notation  used  for  infinite  series  will  be 


S 


t 

n 


(49) 


where  the  n  subscript  denotes  the  n-dependence,  the  k  subscript  denotes 
the  term  and  the  i  superscript  will  be  used  to  distinguish  different 
series.  The  mth  partial  sum  will  be  denoted  by 


k=  1 


and  the  remainder  after  m  terms  as 


R 


l 

nm 


s* 

Jnm 


(50) 


(51) 


4.2  FIRST  SERIES 
4.2.1  Ratio  Test 

Define 


ank  = 

B2k  x2k 

(52) 

(2k  +  n)  (2k)! 

and  for  purposes  of  comparison  define 

al  L  — 

2x2k 

an  k 

(27T> 2k  (2k  +  n)  (]  -  21-  2k) 

(53) 

Q 

Performing  the  ratio  test  on  one  obtains 

an,k  +  1 

r  2k  +  „“i  ri-  2i-2kirxT 

(54) 

an,k 

and 

L2k  +  n  +i|  Ll  -  2— 

7T  _ 

■ — ? 
to  m  W 

"■  + 

l _ 1 

II 

— T 

_2nJ 

(55) 
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Thus,  converges  for  |  x  |  <  2ir.  But  by  Eq.  (48) 

1  ank|  ^  ank 


(56) 


and  thus,  by  the  comparison  test,  Sn  converges  absolutely  for  |  x  |  <  2n. 
This  result  can  be  used  for  obtaining  a  value  for  an  upper  limit  for 

Rnm  I  •  Defining 


(57) 


the  upper  limit  obtained  from  the  ratio  test  is 


R‘r 


2r 


m  +  1 


(2m  +  n  +  2)  (1  -  2-1  “  2m  )  U  -  r) 


(58) 


However,  a  better  upper  limit  can  be  obtained  from  the  alternating 
series  test. 


4.2.2  Alternating  Series  Test 


It  is  noted  from  Eq.  (48)  that  is  an  alternating  series.  The  dif¬ 


ference  between  the  absolute  values  of  sequential  terms  is 

[B2k|  *2k 


I  ank|  -  K.k+l  |  = 


(2k  +  n)  (2k)! 


|B2k  +  2| 

(2k  +  u  +  $)  (2k  +  2)! 


(59) 


Replacing  the  first  term  by  a  smaller  value  and  the  second  term  by  a 
larger  value  obtained  from  Eq.  (48)  results  in 


lankl  —  lan,k+l  I  > 


2rk 


2rk+  1  /  1  \ 

k  +  n  +  2)  \(  1  —  2-1  -2k)!  / 


2k  +  n  (2k 

which  becomes,  after  factorization  of  the  right  side, 

'*w  -  1 » t1  -  ter)  (terr) ']  te 


Noting  that 


and  that 


2k  +  n 


2k  +  n  +  2 


2.2k  +  1 


<  1 


<  -L 


2  2k  +  1  _  1  —  7 

for  all  k,  the  difference  of  Eq.  (61)  will  be  positive  if 


or  certainly  if 


•<f 


X  <  5 


(60) 


(61) 
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Thus,  from  the  alternating  series  test,  an  upper  bound  for  |  R„m  | 
this  range  of  x  is 

I  B  9m  4.  9 

ini  |  j,  I  2m  +  2  ,  A 

l“nm|  <  ; - - — - - 

(2m  +  n  +  2)  (2m  +  2)! 


for 


(62) 


4.2.3  Upper  Bound  for  the  Absolute  Value  of  the  Relative  Error 


The  absolute  value  of  the  relative  error  resulting  from  using  a 
partial  sum  instead  of  the  series  in  Eq.  (46)  is 


IEUI  = 


x"  Hi. 

D‘„  <x) 


(63) 


Any  expression  obtained  from  Eq.  (6  3)  by  replacing  the  numerator  by 
something  greater  and  the  denominator  by  something  smaller  will  be 
an  upper  bound  for  the  absolute  value  of  the  relative  error.  One  obvious 

replacement  is  to  replace  |  Rnm  I  by  its  upper  bound.  Another  replace¬ 
ment  would  be  to  replace  |B2m+2|  by  its  upper  bound  obtained  from 
Eq.  (48).  These  replacements  result  in 


|R‘n, 


2rm  +  1  xn 


Cm  +  n  +  2)  (1  -  2“1 _  2m ) 


(64) 


Since  is  convergent  for  |  x  |  <  2i r,  it  can  be  grouped  in  any  manner. 
Thus,  grouping  as 


Sn  =  (a‘„i  +  a{,2)  +  (a‘n3  +  a^)  +  .  .  .  .  (65) 

and  noting  that  the  first  term  of  each  group  is  positive  and  the  second 
term  is  negative,  from  Eq.  (61)  it  is  deduced  that  each  group  is  posi¬ 
tive  for  x  <  5.  Thus,  is  positive  and  one  obtains 


D*„(x)  >  xn  Pi-  _  —L__] 

|_n  2(n  +  1)  J 

Using  Eqs.  (64)  and  (65)  to  replace  the  numerator  and  denominator, 
respectively,  of  Eq.  (6  3)  one  obtains 


But 


El 

nm 


< 


2  2m  +  2  j.m+1 


(22m  1  _ j)  (2m-„  +  2) 

1  X 

_n  2(nH  1) 

(67) 


(2m  +  n  +  2)  P  J_  _  x  1 

Ln  2(n  + 1)  J 


(2m+n  +  2)  (2n  —  nx  +  2) 
2n(n  +  1) 
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and 

and 


so 


2n~nit  +  2  =  2-x  +  2-  >  2-x 

n  n 


2m  +n  +  2 
n  + 1 


>  1 


IE*. 


< 


22m+3 

22m+l_j 


fm+l 

2-x 


(68) 


This  upper  limit  approaches  zero  as  m  approaches  infinity  and  is  easily 
calculated.  Note  that  it  is  also  independent  of  n  and  thus  holds  for  all 
orders. 


4.3  SECOND  SERIES 


4.3.1  Integral  Te»t 


The  terms  of  the  second  series  are 


a’  =  2 


nk 


j=0  (n-j)  ! 


xn-je-k* 


(69) 


It  is  shown  that  a^k  satisfies  all  of  the  conditions  for  the  integral  test. 


o 

Thus  an  upper  limit  for  Rnm  can  be  obtained  as 


ir™'  <1^-  . 


n-i  /  e“xt  rj_1  dt 


(70) 


Integrating  successively  by  parts,  one  obtains  for  j  >  0 

J*  e  “xt  t“i_1  dt  =  [e_xtJ  £(-l)  y  "  -  [«.-][«»  j  (j-l)  J 

m  m 

+  [<-x>’«-]  [<-»■  7^5r] 

+  (-l)i_1  ^(-x)j-1  e"xt]  [(-1)1  -Llj 

OO 

(-l)j  (-x)j  (~I)j  f  e~xt  dt 

i  I  J  t 


21 


AE  DC-T  R-68-254 


or  after  simplification 

30 

J*  e~ *l  t“i+1  dt  =  2  (-I)*"1  ^7^ !  -iULfl 


i-i 


i=l 


e— rax  +  -  E(mx)  (71) 


where  E(mx)  is  the  exponential  integral  discussed  in  Ref.  5.  Substi¬ 
tuting  Eq.  (71)  into  Eq.  (70)  results  in 

lR‘- 1- i  L  b"1  gffir 

+  f=1  (-1)J  j  i  (n_j)  i  x“  E(mx^>  +  TT  X“  E(mx)  (72) 

Note  that  the  last  term  can  be  included  in  the  last  summation  by  letting 
the  range  of  the  index  j  start  at  zero.  Making  this  change  one  finds 

£^(-1)*  -7^n_  ■),  xnE(mx)  =  xnE(mx)  21 

The  resulting  summation  is  zero  as  can  be  found  in  Ref.  5.  Thus,  the 
right-hand  side  of  Eq.  (72)  reduces  to  the  double  summation.  This  can 
be  simplified  to 


|R 


nm| 


<  i  i  i-i),"(")(i  -  l)!i7r 


(73) 


i  —  1  j  =  i 


It  can  be  proved  that  the  double  summation  of  Eq.  (73)  is  just  a  re¬ 
arrangement  of  the  terms  of  the  double  summation  of  Eq.  (72).  Simpli¬ 
fying  further 


IRnml  <  e- 


But,  from  Ref.  5 


MJ'U-mxSji  £(-!)■  ft)  (74) 

1—1  J  =  i  \/ 

(“)  -  C"=0  *Cl‘) 0  <  i  <  “ 


so 


W)1  (■)  ■  I  «*(;:/)*  l/-1'1  (7  V  <-»■ 

=  21  <-Di+1  ("71)  +  L  (-Dj  (n“l)  +  (-i)n 

j=i-l  j  =  i  ^  ' 

=  n£'[(-i)j+i  +  (-Di]  (n-')+  (-d1^-;)  +  (-D"-1  ("i;)  +  (-Dn 
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It  is  readily  shown  that  the  summation  is  zero  and  that  the  last  two  terms 
cancel  leaving  only 

(75) 

Substituting  Eq.  (75)  into  Eq.  (74)  one  obtains 

n 

IR’nml  <  e"mJt  I  (i  - 

i=  1 

or  finally 

|RJnml  <  (n  -  1)  !  ®— * 


1) 


,  (D~l\  Jilli 

V  “  1/  mi 


n  —  1 


y  x» 

i  =  0  i!mn-i 


(76) 


4.3.2  Upper  Bound  to  the  Absolute  Value  of  the  Relative  Error 

The  absolute  value  of  the  relative  error  caused  by  using  a  partial 
sum  instead  of  a  series  in  Eq.  (47)  is 


I  Enml  = 


Da„  (*) 


(77) 


Any  expression  obtained  from  Eq.  (77)  by,  replacing  the  numerator  by 
something  greater  and  the  denominator  by  something  smaller  will  be  an 
upper  bound  for  the  absolute  value  of  the  relative  error.  The  obvious 
replacement  for  the  numerator  is  the  upper  bound  just  obtained  for 
2  I 

Rnm  I .  This  bound  will  be  further  increased  in  value  by  replacing 


nm 

each  term  in  the  summation  of  Eq.  (76)  by  a  greater  value,  specifically 

n-  1 

l 

i  =  0 


|R?n 


<  (n  -  1)!  e“mj 


x1 

i! 


(78) 


Since  all  of  the  terms  of  the  double  summation  of  Eq.  (47)  are  positive, 

n 

then  Dn(x)  is  greater  than  the  first  n  terms.  Thus  after  reindexing  it 
is  found  that 


n  —  1 


Dn (x)  >  £  y  e"*xj 


i=  0 


thus,  obviously 


n  —  1 


DJn(x)  >  (n  -  1)!  e“x  £  ~ 


i  =  o 


(79) 
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Using  Eqs.  (78)  and  (79)  to  replace  the  numerator  and  denominator, 
respectively,  of  Eq.  (77),  one  obtains 

|EU  <  (80) 

This  upper  limit  is  easily  calculated  and  approaches  zero  as  m  approaches 
infinity.  Again  the  bound  is  independent  of  n  and  thus  holds  for  all  orders. 


4.4  RANGES  OF  CALCULATION 

The  upper  bound  for  the  absolute  value  of  the  relative  error  is  re¬ 
lated  to  the  number  of  significant  figures  by  a  theorem  proved  in  Ref.  8s 

If  the  absolute  value  of  the  relative  error  of  any  number  is  less  than 
5  x  then  the  number  is  certainly  correct  to  s  significant  figures. 

Using  this  relation  the  accuracy  of  each  approximation  was  plotted 
and  is  shown  in  Fig.  4.  Given  a  value  of  x  and  the  number  of  significant 
figures  desired,  from  Fig.  4  one  can  find  which  approximation  and  what 
order  is  required.  Thus  one  can  construct  a  table  dividing  the  calcula¬ 
tion  into  ranges  of  x  and  the  best  approximation  to  use.  Table  IV,  for 
example,  gives  such  a  division  for  five  significant  figures. 

It  is  shown  that  the  approximation  for  each  range  falls  below  the 
line  s  =  5  in  Fig.  4.  Since  the  small -x  approximation  of  order  5  and 
the  large-x  approximation  of  order  8  intersect  below  this  line,  one 
never  needs  to  go  higher  than  these  orders  to  obtain  five  significant 
figures.  Ranges  13  and  14  were  determined  by  the  limitation  of  the  IBM 
System/ 360  computer  that  the  order  of  the  exponential  function  be  less 
than  174. 

A  subroutine  was  written  to  calculate  the  Debye  function  to  five  sig¬ 
nificant  figures  using  the  ranges  of  Table  IV.  The  subroutine  is  called 
by  the  statement: 

CALL  DEBYE  (N,  X,  F,  G) 

The  arguments  N  and  X  input  n  and  x,  respectively,  where  the  order  of 
the  Debye  function  must  be  in  the  range  0  <  n  <  6.  The  Debye  function 
and  its  complement  are  returned  by  the  arguments  F  and  G,  respec¬ 
tively.  The  listing  of  this  subroutine  is  given  in  Appendix  II. 


24 


Fig.  4  Accuracy  of  Approximations 
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TABLE  IV 

RANGES  OF  CALCULATION 


Range  No. 

Range 

a. 

m 

1 

0  <  x  <  0.  007 

l 

0 

2 

0.  007  <  x  <  0.  22 

l 

1 

3 

0.  22  <  x  <  0.  68 

l 

2 

4 

0.  68  <  x  <  1.  12 

l 

3 

5 

1.  12  <  x  <  1.  52 

l 

4 

6 

1.  52  <  x  <  1.  80 

i 

5 

7 

1.  80  <  x  <  2.  03 

2 

8 

8 

2.  03  <  x  <  2.  44 

2 

7 

9 

2.  44  <  x  <  3.  05 

2 

6 

10 

3.  05  <  x  <  4.  07 

2 

5 

11 

4.  07  <  x  <  6.  10 

2 

4 

12 

6.  10  <  x  <  12.  21 

2 

3 

13 

12.  21  <  x  <  87 

2 

2 

14 

87  <  x  <  174 

2 

1 
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APPENDIX  I 

FORTRAN  LISTING  OF  MAIN  PROGRAM 


C001 


0002 

0003^ 

0004 

"0005 
0006 
0007" 
0008 
0009 
0010 
0011 
0012 
00  l-j" 
0014 
0015 
_  C016 
0017 

_ 0018 

0tfl9 

0020 

0021" 

0022 

0023' 

0024 


0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

O'04O 

0041 

0042 

0043 


C  AEA00052  BL ACKBGDY  EMISSIUN  D.C.  TODD"  2-12-68  SY8CC2-VOO' 

REAL  UNIT!  5*41/' H  ICR* i*  CH  H  IN  •_,*  FT  S  ««»  HINS 

l"  "  •  HR  2*0., ■KDEG*  ROEG' i'CDEG' t  •  FOEfc*  ,£).,•  ERG*.*  J  ', 

2  '  C AL • » • FTLB • « 1  BTU* /, D 1ST! 5) /I ,E-6, .0  It  1. • .0254, . 3048/, 

....  ,  7j  ",60 .7 3600 , TEflPt 4 ) /l  .  ,  .  555555671 .  , . 5555556/  , 

4  TURE 141/0. ,0. »273. 15,459.67/ ,ERGVl 5) /l.E-7,1. ,4. 19002,1.35582, 
"5  ”  1  O' 5 5.8  7  /  ,'C  R  /'2  .  9  97925  E8  /  ,  HR/'6\'62 56E'-  3  4  /  ,  K  R71718U5  4  E-2  3 / , 

6  P 1/3. 14 1593/, ZE TA3/ 1.20 2057/, ZETA4/1. 082323/, X2/ 1.593624/, 

7"  X 3/ 2.82 1439/ «  X4/ 3 ."92069/^X5/4. 965 114/*  " 

8  K,LAH,NU,NHE,M1F,LPE,LHF,IC0H 

- - 

WRI TE 16 i 1000 ) 

TREAD!  5,TO"Or."END=“9CdViW,  (L,lS,~lT;TE,LTfitT,ao,DQ 
C  UNIT  CONVERSIONS  AND  CALCULATI_ON_OF^  PHYSICAL  CONSTANTS_ 

.  IF(  IW.LE.O  I  GO  TO  41"  “  . 

A=DI  ST  I  1HI/0ISTI  ILI 
■  C-CR4  TTM'£‘(  IS  IT'D  I  ST  fi n 

H=HR/ IERGY I IEI *T  IHEI  IS )  I _ 

. KVKR*TE'H'P'(TTI/ERGY(I  EJ 

TREF=TURE< IT! 

'  KO'H=K/H  . . 

TAU=4.*PI*(KOH/C)**2*KOH 
'  S I  G=  3 .  *Z  E  T  A  4  *K  *TA  U 

_iiyy= AU _ _ 

B2=K0hVX2  .  v 

B3=K0H*X3 

B4=C/I  A*K0H*'X41 . . 

B5=C/(A*KOF*X5» 

Cl=2".*T>[*H*ir*C/A**4 
C2-C/ 1  A*K0H) 

'  "  IF  (  I  P.EQi'  1 JGO"  TO  21 

I  P=  1 

"  '  HR  (TE  (  6 ,1005  I 

21  WRITE(6,1002)UNIT(IW.1),UNIT(IL,1>,UNIT( IS, 2), UNIT! IT,3I, 

- I - uNIttt  E,4)  ,A ,C,H,K,TRFF - 

HRITE(6,1003ISIG,TAU,B2,B3,B4,B5,C1,C2 
C  CALCULATION  "OF  TEMPERATURE  DEPENDENT  QUANTTTTE3 

41  TCK=T  +TREF 

'  '  I'FITCK.LE'VO.IGO  TO  61 
TA=TCIC 

E  B=  S I  G*  T  A*"T  A  *  fA*T  A 
FB=TAU*TA*TA*T  A 
NME  =  B3*T"A 
NHF=B2*TA 

' . .  LME=B5'/TA"  ' 

LMF=84/T  A 

- T  (Terr.  E'0'."2)  Gd"T0_4"2" 

...  I.P.=  2_  _ _ .  .. _ _  _ 

MR  IT  E 1 6 , 1006  I 

42  WRITE(6,1002)T,TA  ,EB  ,FB,  NME  .NHF  ,  LME  >LMF 

'  C  CALC ULAT ION  OF' "WAVELENGTH  DEPENDENT  QUANTTlTK 
61  1 F ( N*  LE. 0) GO  TO  1 

- . gotot72;7  3',74t;l' . 

72  0N=C/(150.*A*KGH*TA) 

.  UM="C/T."00"1*"A*KGH'*TA) 

GO  TO  75 
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APPENDIX  I  (Concluded) 


0044 

73 

dN=.O0L*KOF*TA 

004b 

QH=15Q.*K0H*TA 

0046 

GO  TC  75 

C047 

74 

UN=.001 

CG4  8 

UM=150. 

0049 

75 

CONTINUE 

OObO 

DO  66  1=1, N 

0051 

0=00+11-1)  »ca 

00b2 

IFlU.LT.CN) c=cn 

0053 

IF(C.GT.QM)C=CF 

CC54 

GO  TO!  62,6  3, <4  1,  L 

0055 

62 

LAH=C 

0056 

NU=C7 1 A*LAH ) 

0057 

X=.JU/IKOH*TA) 

0058 

GC  TC  65 

C059 

63 

NU  =  Q 

0060 

X=NU/1K0H*TA> 

0061 

LAM=C/(A*NU) 

0062 

GJ  TO  65 

0063 

64 

X=3 

0064 

NU=KOH*TA*X 

0065 

LAW=C/1A*NL) 

0066 

65 

FN=2.*PI*(NU/C)**2/IEXP(X)-1.) 

0067 

EN=H*NU*FN 

0068 

EL=NU*EN/LAM 

00o9 

fl=nu*f\/lam 

00  70 

CALL  CFBYE(3,X,EBN,E8L) 

0071 

EdN=EU*CBN/16.*ZETA4> 

0C72 

£BL=E8*EBL/16.*2ETA4) 

00/3 

CALL  CEBYEI2,X,FBN,FBL) 

00  74 

FBN=FB*FBN/(2.»2ETA3) 

00  75 

FBL=F8*FBL/12.*ZETAJ| 

00  76 

IF  IIP. EO. 3) GO  TO  66 

0077 

I  P  =  3 

CO  78 

WRITE (6, 10C7) 

CC  79 

66 

WRITE  16, 100  •  )LAM,NU,X,EN,FN,EL,FL,EBN,FBN,EBL,FBL 

008C 

GO  TO  1 

OOd  1 

900 

HRITEI6,  10C41 

0082 

STOP 

CCJ3 

1000 

FGKWAT (• 1AEAC0C52  O.C.  TUDD  BLACXBODY  EMISSION1! 

0084 

1001 

FORMAT (6 11 ,I6,JE12.3> 

C0d5 

1C02 

FORMAT  1 LH03XA4 ,4 1 3XA4] ,4X  LP5E12.4] 

0086 

1003 

FOrtMATl IP11F12.41 

C03  7 

1004 

FORMAT!  •  CTFE  ENO'/IHL) 

CO  18 

1005 

FJRMAT  (  1H04X2HH1CX2HIL10X2HI  S  10X2H I T  10X2HI  E  il'X  1HA  riVfHCTT'XlHH  1  lx 

1 

IHK9X4HTMEF/4X5HSI'J‘'A3X3HTAU9X2HB210X2HB310X2HB410X2HB5 

2 

!  10XZFC l 10X2RC2I  . 

008  9 

1006 

FORMAT! 1H05X  1HT 1 0X2H TA 1C X2HEB 1 0X2HF B10X3HNME9X3HNMF9X3HLME9X3HLMF 

1 

/ )  ~  - 

0040 

1007 

FOR«AT 1 1H02X6HL  AME0A8X2HNUI IX1HX10X 2HEN1 0X2HFN1 C X2HEL 10X2HFL 10 X 

1 

3HF3N9X  3HFBN  5  X3HE  8L  9  X3HFBL / ) 

0081 

END 
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FORTRAN  IV 


0031 

0002 


000  3 

0004 

CCC5 

0006 

0007 

OOOS 

0009 

0010 

COLI 

00L2 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
_002l 
0022 
0023 
0024 
0025 
0026 
0027 
C  02  8 
0029 
0030 
0031 
0032 
0033 
0034' 
0035 
0036 
0037 
0038 
0039 
0040 


APPENDIX  II 

FORTRAN  LISTING  OF  SUBROUTINE  DEBYE 


G  LEVEL  1>  MOO  2  DEBYE  DATE  ='68213 .  19/32/19 

C  AEA00048  DEBYE  FUNCTION  O.C.  TODD  1-18-68  ST8OC2-V0O 
SUBROUTINE  CEBYE I N, A , F , G I 

REAL  B2(  5)/.  1666  667, -.03333333,  .C2380952,-.  03333333  ,.  07575758/ 

1  ZET All  6)/ 1.6449 34,  1.202057,1.08232  3,1.036928,1.01734  3,1.008349/, 

2  FAC  I  10)/l.,2.,6.,24.t  12C.  ,  720.  ,504  0.  ,40  320  .  ,362  8  8C.V362  88750:7; 

3  A(5  )/. 007, .22, .68, 1.12, 1.52/, 

4  BI81/174., 87. ,12. 21, 6.1,4.C7,3.C5, 2.44,2.03/'  "  . 

IF (X.GT.l.S)GO  TO  20 

C  F  CALCULATED 
M=0 

'  DO  11  J- 1,5  " 

IFIX.LE.AI J) IGC  TC  12 

11  M=M+1 

12  XS=X*X 

SUM=l./N-.5*X/IN*H 
IFIM.EQ.OIGC  TO  14 
POW=l. 

KZ=0 

DO  13  K  =  1,M 

K2=K2*2 

POh=POW*XS 

13  SUM=SUH*B2(K)*P0W/I I K2+N I *F AC [ K2 1 1 

"  l4'F=SUM*X*'*N  . .  . . 

G=FACIN)*ZETA1(N)-F 

RETURN 

C  G  CALCULATED 

20  M=0 

DO  21^J=1,B 
IFlX.G^it'j  flGO'  TO  22 

21  M=M+1 

22  IFIM . Gt .0  IGC  TC  2  3 

G=0. 

GO  TO'  26 

23  SUM=0. 

PbR=OV 

DO  25  K=  1 ,  M 
POh=POW*X 

_ SM=1./K**(N+1I  _ 

PW=SM . 

DO  24  J=1,N  . 

.  P'M=PW*X*K .  '  _ 

24  SM=SM+PW/FACIJI 

25  SUM=SUM*SN*EXP<-PCtf] 

G=FAC(N)*SUH 

26' F=FACI N! *ZE T Al IN l-G  ' 

RETURN 

. . END . 
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